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The Construction and Use of Divergence Free Vector 
Expansions for Incompressible Fluid Flow Calculations 


Neatan Mae Giolla Mhuirts 

Institute for Computer Applications in Science and Engineering, 
Mail Stop 132C, NASA Langley Research Center, 
Hampton, Virginia 23665, USA. 


ABSTRACT 

For incompressible fluids the law of mass conservation reduces to a con- 
straint on the velocity vector, namely that it be divergence free. This constraint 
has long been a source of great difficulty to the numericist seeking to discretize the 
Navier-Stokes and Euler equations. In this paper we will discuss a spectral 
method which overcomes this difficulty and we will demonstrate its efficacy on 
some simple problems. The velocity is approximated by a finite sum of divergence 
free vectors, each of which satisfies the same boundary conditions as the velocity. 
Projecting the governing equations onto the space of inviscid vector fields elim- 
inates the pressure term and produces a set of ordinary differential equations that 
must be solved for the coefficents in the velocity sum. The pressure can then be 
recovered if it is needed. 


The leaearch for the author was supported in part by NSF Grant 830713 and by the National Aeronautics and 
Space Administration under NASA Contract No. NASI - 17070 while he was in residence at ICASE, NASA 
Langley Research Center, Hampton, VA 23665-5225. 
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1. Introduction. 


We will be interested in the flow of incompressible fluids for which the physical law of mass con- 
servation reduces to a constraint on the velocity vector of the fluid, namely that it be divergence 
free. The pressure in such flows is not then a thermodynamic variable determined by an equation 
of state but rather can be thought of as a Lagrange multiplier which adjusts itself instantaneously 
to ensure that this kinematical constraint on the velocity vector is satisfled. There is no evolu- 
tion equation for the pressure nor does it satisfy any predetermined boundary or initial condi- 
tions. 

Numericists, seeking to solve the governing equations approximately, have found that their 
greatest difficulty lies in the treatment of the pressure variable. While many ingenious methods 
have been devised to overcome the difficulties, the treatment advocated in this work is in a 
mathematical sense the most natural and offers many computational advantages. Here, the pres- 
sure term is eliminated from the equations entirely, and the divergence free condition is satisfied 
exactly by the numerically obtained approximation to the velocity vector. Moreover, as the com- 
ponents of the velocity are expanded in terms of series of polynomials which are solutions of a 
singular Sturm Liouville problem, whose excellent approximation properties are well documented, 
(e.g., Gottlieb & Orszag [1977], Quarteroni [1983]) convergence of our approximation will be 
bound only by the smoothness of the solution and by the number of terms used in the com- 
ponent expansions. For infinitely differentiable velocity fields we should expect to achieve 
"exponential convergence" (Canute et al.) . 

The essence of the method involves expanding the velocity in a series of divergence free vec- 
tor fields each of which satisfies the same boundary conditions as the velocity and where in gen- 
eral the coefficents in the expansion are time dependent. The infinite sums are truncated and 
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substituted into the governing equations, which are the Navier-Stokes or Euler equations. Inner 
products are taken with vectors fields that satisfy inviscid boundary conditions, eliminating the 
pressure term and giving rise to a set of coupled ordinary differential equations for the coefTicents. 
Nonlinear terms are treated explicitly in the time marching scheme, but the linear viscous terms 
can be treated implicitly. The method can also be used to solve linear stability problems; for 
these we will get a generalized matrix eigenvalue problem rather than a set of ordinary 
differential equations. The eigenvalue determines the stability of the base flow and each eigen- 
function is the set of coefficents in the expansion of the corresponding perturbation velocity field. 

In the next section of this paper we will see how to construct complete sets of solenoidal vec- 
tor fields assuming that these are periodic in two directions. There have been various applica- 
tions of the method both to simulate certain fluid flows (e.g., Moser & Moin’s [1984] work on the 
infinite Taylor-Couette system) and also to solve linear stability problems for rather complex 
base flows (Mac Giolla Mhuiris [1986]). For pedagogical purposes, in this short paper we will con- 
sider the linear stability of a rotating Poiseuille flow. Extensive applications of the method to a 
set of more taxing stability problems that arise in the context of a fluid flow phenomenon known 
as vortex breakdown will be published elsewhere. We will also indicate how the method was used 
to simulate the infinite Taylor-Couette system. 


2. Divergence Free Vector Bases. 

It is well known (Ladyshenskaya [1966]) that the space of square integrable vector functions, 
Zr^(D) defined on a bounded domain D [D C jR"’ n = 2,3) can be decomposed into those that are 
divergence free and whose normal components vanish on the boundary and those that can be 
expressed as the gradient of a differentiable function defined on D. The domains of interest for 
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this paper are either cylinderical or annular. For now we will consider vector fields defined on the 
section of a cylinder, T which are periodic in both the axial and azimuthal variables, having as 
their axial period the tube length, 2H. 

We will decompose L"^[T) as follows. 

L^{T) = J(T) + /(T), (2.1) 

where. 


J{T] = 


(a) V'j£ = 0 T, 
u e L^{T) (6) u’n = 0 on dT, 

(c) «| = «! 5 , • 


( 2 . 2 ) 


and ^2 represent the ends of the cylinder. Given in this form J{T) is clearly the space of 
(a) incompressible, (b) Inviscid, (c) periodic velocity fields. 

The set of "viscous" velocity fields on T is a subset of J{T) denoted J^{T). 


J^[T) =\ueJ{T) I u = 0ondT 


(2.3) 


An alternative representation of J{T) (Richtmyer [1978]) is given by, 

(a) <u, VP> ~ 0 /or P f C^[T) 


J{T) = 


U€L\T) 


(6) u| = u| 5^ 


where T is the closure of T and <•,•> represents the usual inner product in L^{T), 


(2.4) 


<u,v> = u'v rdrdddz 


(2.5) 


The space J{T) endowed with this inner product is a Hilbert space and a closed subspace of 
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L^T). 

The projection of L^{T) onto J{T) will be denoted by IT. It is clear that vectors of the form yp 
are perpendicular to all u in J(T) and in fact (Ladyshenskaya [1966]), 


J^[T) =1 ueL^{T) I ^ = '^p for some p in C^{T)j . (2.6) 

n then has the following properties, 

Il:L^{T)-*J{T), ( 2 . 7 ) 

Il£=u for all u e J(T), (2.8) 

n yp = 0 for all p e C^(T). (2.9) 


To solve the Navier-Stokes equations in the domain T with periodic boundary conditions in 
the axial and azimuthal directions, we seek a vector u in J°(T) satisfying, 

^ + (u-y)u = -yp + (2.10) 

where i/ is the kinematic viscosity of the fluid. 

We can eliminate the yp term in this equation by means of the projection IT and, using 
some vector identities, we arrive at the following equation for u. 

^ = n| wxu - i/yxwj (2.11) 

where w is the vorticity, w = yxu. 

Projection of the Navier-Stokes equation onto a finite dimensional subspace, J^[T) of J{T) 
is achieved in practice by taking the inner product of the equation with basis vectors for 
This process eliminates the operator II from the equation, for given any vector / in L^[T) and 
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any vector A in J{T) we have that, 


<n/, A> = </, A>, (2.12) 

as projections are self adjoint and as the projection of any vector in J{T) is itself. 

It is worth emphasizing that although we are solving the equations for a viscous fluid we 
nevertheless project the governing equations onto the space of inviscid vector fields. The reason 
for this is that having found a velocity u such that the vector / defined as, 


/ = Uj — wxu + (2-13) 

is orthogonal to all ^in J{T) then / e J^{T) and so there exists a scalar function p (a pressure) 

with / = Vp. If, however, / were in (T), which contains J“(T), the existence of a pressure is 
not guaranteed; consequently, u may not correspond to a physical solution. 

Leonard and Wray [1982] demonstrated a divergence free vector function expansion on a 
cylindrical domain for viscous velocity fields that are Fourier decomposable in both the axial and 
azimuthal variables. Mac Giolla Mhuiris [1986] considered expansions in terms of somewhat 
different basis functions, and it is these which will be described here. 

The velocity field, u, satisfies the continuity equation and is Fourier decomposable in z and 
0, which means in effect that only two of its three components, u, v, w are independent. This 
motivates the introduction of two vector families, xt expansion of the form. 


a = Ei «2n*m(0x;T('-) + a2n-litm(0x„(r) ^ e . (2.14) 

nkm ( I 

The components of the vectors are found as follows. Expand two of the velocity components, 
say the first and the third, independently as. 
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u 




(2,15) 


nkm 


(2.16) 


where /*(r) are complete sets of polynomials chosen to satisfy the boundary conditions that are 
imposed on u, w. The r and z components of x* have now been picked, and it remains for us to 
choose the 0 components in a manner that ensures the vectors xt ^ ^ divergence free. 

Consider for example, Xn(j'). 


V- Xn{r) e = 0, 


=> 


(2.17) 


(»•/«('•))' + i^Xn,e = 0, (2.18) 

where the prime denotes a derivative with respect to r. This equation gives us the 6 component of 
X^(r). Rescaling, it is found that an expansion of the form (2.14) is possible for non-zero azimu- 
thal wavenumbers where, 

X"(r) = ( tm/rTCOi - (’•/,T(0)'> o) . (2.19) 

Xn('-) = ( 0, - rA:/+(r), m/ + (r) j (2.20) 

and such an expansion will guarantee that u is divergence free. This expansion is clearly incom- 
plete for azimuthal wave number zero, (m = 0), i.e., for axisymmetric flows and for that case 
the following expansion vectors can be used. 


X„-(r) 


' 

T 


3 


( 2 . 21 ) 
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= (o, /„+(r),o] . (2.22) 

The polynomials /^(r) must be chosen so that the vector u given by (2.14) satisfies 
appropriate (viscous or inviscid) boundary conditions on the walls of the domain, T, and is single 
valued at the origin, r = 0. We will denote the polynomials used in the inviscid case by a*(r) 
reserving /^(r) for viscous expansions. We have then upon truncating (2.14) an approximation 
to u of the form, 

N K M 

= XI Ij Ij “nltm(0 

n=l k=—Km=—M 

where, 

^ntjr,z,0) = ^(r;A:,m)e‘(*^ + "‘'). (2.24) 

The projection vectors will have the same form as the expansion vectors, i.e., we will project 
with vectors, ,0) , where 


for 




= ef (r; k, m) e 


(2.25) 


/ = 1,...,W; p = q = 


with the vectors ^ being given by equations (2.19 - 22) using the inviscid polynomials, af'{r) for 
the components. 


It is possible to choose the polynomials and /*(r) in many different ways. Leonard & 

Wray [1982] in their consideration of certain turbulence simulations employed an unusual class of 
Jacobi polynomials to reduce the bandwidth of the final matrix system. These polynomials were 
also used by Spalart [1983] in his simulations of boundary-layer transition. Moser &; Moin [1984], 
in their work on the infinite Taylor-Couette system, used Tchebychev polynomials and incor- 
porated the weight function, against which these polynomials are orthogonal, into the projection 
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vectors. Mac Giolla Mhuiris [1986], investigating the temporal stability of a set of flows that are 
said to model vortex breakdown, (Leibovich [1984]) experimented with various sets of polynomi- 
als and found that Legendre polynomials were, perhaps, the easiest to use. All of the above sets 
are solutions to singular Sturm Liouville problems, and consequently we can expect expansions in 
terms of any of these polynomials to exhibit excellent convergence properties. 

The single valuedness criterion, which must be applied along the center line of the tube for 
the vector u, causes the polynomials /*(r) and a/^(r) to depend on m, the azimuthal 
wavenumber (Joseph [1970]). One appropriate choice for a*(r) in the case of pipe flow is. 


a,+(r) = rP,(2r-l) 


for all m, 


aj (r) = (1 - r)P,(2r - 1) if | m| = 1, 


a-(r) = r(l - r)P,(2r - 1) if I "^1 


where the radial variable has been scaled by the tube radius and P/(r) is the Legendre polynomial 
of order I (Abramowitz & Stegun [1970]). The corresponding choice for f^{r) is. 


/+(r) = r(l - r)P„(2r - 1) 

/n (0 = (1 - rfPni'^r - 1) if 1 m| = 1, (2.27) 

/-(r) = r(l - r)2p„(2r - 1) if I »n| 1. 


For the problems tackled to date, all of which separate in some sense in the azimuthal direc- 
tion, this dependence on m presents no difficulty. We solve separate problems for each azimuthal 
wavenumber; so having chosen an m the expansion and projection sets are fixed throughout the 
calculation . Indeed, in theory there is no difficulty even if the problem at hand is truly three 
dimensional; however, some care is required in implementing the method to ensure that the 
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correct radial polynomial set is being used for each azimuthal component of the velocity. 


3. Linear Stability Problems. 

To determine the linear stability of a solution U of, say, the Navier-Stokes equations, we 
consider whether infinitesimal perturbations to that solution grow in time. Therefore we linearize 


the governing equations about U and seek solutions of the form, 

u(r,^,0)e--‘. (3.1) 

The character of a then determines the linear temporal stability oi U. If <r = a + where a, ^ 
are real, then 


/3> 0 

= > 

U 

is unstable 

/?=0 

= > 

u 

is neutrally stable 

/3< 0 

= > 

u 

is stable. 


The equations that must be solved have the form, 

iVu = Eu + Re 

E and S are operators defined on J{T) as follows, 

Su = -n(vx5i^ 


(3.2) 


(3.3) 


(3.5) 


and 

Eu = Yl{t£xU + nxu), (3.6) 

where Cl is the base fiow vorticity, Q = yx^. Some suitable nondimensionalization has intro- 
duced the Reynolds number. 
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Re= 



(3.4) 


Rq and Uq being characteristic length and velocity scales associated with the base flow. We can 
take Rq to be the radius of the tube and Uq to be the maximum value of the base flow azimuthal 
velocity. 


We will solve (3.3) approximately by using the expansion for u, taking inner products 

with the projection vectors, and solving the resulting generalized matrix eigen problem for 
the eigenvalues <r and the eigenvectors (In this case, of course, these coefficents do not 

depend on time.) In general this matrix problem can be written as, 

(3.7) 

The matrix A is purely real and arises from the fact that the expansion and projection vectors 
are not orthogonal. 

^Inpkqm = <^pq>Rnkm>i = < Xn>^pk^qmJ = ^\n^pk^qmi (3.8) 

The Kronecker delta symbol, arises because the Fourier bases employed in the axial and 

azimuthal directions are orthogonal. The matrix B arising from the convection terms is also 
purely real. 

^Inpkqm ~ "^^pqt (V^:^itm)^^ + ^l^^nkm^' (3.9) 

Finally, the matrix C arising from the viscous terms is purely imaginary. 

^Inpkqm ~ ^ ^pq’ • (3.10) 

Using the orthogonality of the Fourier bases, it can be written as. 




- 11 - 


^Inpkqm ~ ^\n^pk^qm. (^-H) 

The form of the matrix B depends on the base flow U. If U is independent of 2 and 0 then it will 
be possible to find a matrix B such that, 

(3.12) 

This simply corresponds to the fact that the partial differential equation (3.3) separates in z and 
6, for base flows that are independent of those variables. The stability of such flows can be deter- 
mined by solving the NxN generalized matrix eigen problem, 

aAa = \ B + -^C a. (3.13) 

Re 
V 

The matrices A, B, and C all depend parametrically on the wavenumbers of the projection and 
expansion vectors so it is advisable to separate them into submatrices that can be evaluated 
independently of these and any other parameters occurring in the base flow. The submatrices are 
evaluated once and then stored in the computer. The required integrations can be done at very 
little cost by utilizing the orthogonality properties of the expansion and projection polynomials. 
The full matrices can then be reassembled without the need for doing any further integrations. 

One can always band the A and C matrices by appropriate choice of the polynomials /„ (r) 
and a*(r). However the matrix B will generally be full, though for certain rather simple base 
flows such as the Hagen Pouiseille flow considered be Leonard & Wray [1982] it is also possible to 
band B. It is also worth noting that the matrix problem we get when considering the inviscid 
stability of base flows is a purely real one, and consequently the eigenvalues occur, as they should 
do, in conjugate pairs. 

To demonstrate the method let us consider the viscous stability of a rotating Poiseuille flow. 
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E = ( 0, Vgr, (1 - r^) J . (3.14) 

Cotton et al. [1980] found that this flow with Vq = 0.2147 was neutrally stable to disturbances 
having azimuthal wavenumber, m = 1 and axial wavenumber, Jt = -1. The following table lists 
the most unstable eigenvalue found for this flow using a code which implements the method we 
have discussed in the previous section (with Legendre polynomials as the radial basis functions). 
The first column of the table gives N, the number of radial basis vectors that were used to obtain 
the eigenvalue given in the next two columns; N is also the order of the matrix problem that 
needs to be solved at each step. 


Most unstable eigenvalue found for the rotating Poiseuille flow (3.14) 
with m = 1, k = -1, and = 0.2147. 

N 

frequency 

growth rate 

4 

-0.00029 

.00334 

6 

-0.00279 

.00101 

10 

-0.00284 

.000001 

14 

-0.002847 

.0000001 

18 

-0.002847898 

.0000001379 

22 

-0.002847898 

.0000001378 


The convergence is exponential in N or some power of iV, and there is no evidence of significant 
roundoff error. 

We also used the divergence free expansion method to solve the adjoint stability problem 
which can be written as follows, 

iXu = E*u 4- (3.15) 


The operator E is the adjoint operator to E and is given by. 
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E*u = -n| Uxu + Vx(h^^) • (3.16) 

The direct and adjoint spectra obtained by solving (3.3) and (3.15) should, of course, be conju- 
gate to each other, and how well a numerical scheme reproduces this theoretical result is a test of 
its accuracy. The following table lists the corresponding adjoint eigenvalue found by solving 
(3.15) for the base flow and wavenumber pair considered previously. 


Eigenvalue found by doing the adjoint problem (3.15) for the flow (3.13), 
with m = 1, k = —1, and Vq = 0.2147. 

N 

frequency 

growth rate 

4 

-0.00270 

-.000094 

6 

-0.00280 

-.000013 

10 

-0.00284 

-.000004 

14 

-0.002847 

-.0000002 

18 

-0.002847898 

-.0000001379 

22 

-0.002847898 

-.0000001379 


Clearly the agreement between the adjoint and direct results is excellent. 


4. Flow Simulations 

The global nature of the expansion functions we are considering here causes them to be extremely 
inefficient for the direct evaluation of nonlinearities such as the convective derivative term that 
occurs in the Navier-Stokes and Euler equations. However it is possible to use the so-called collo- 
cation technique (see for example Orszag [1972] or Kreiss & Oliger [1971]) in conjunction with 
our divergence free expansions to alleviate this problem. 

The idea is to evaluate the nonlinear terms in physical space, expressing the result of this 
calculation (a vector function known on a set of grid points), componentwise, as a series of the 










-14- 


selected basis polynomials being used to form the projection vectors. The orthogonality proper- 
ties of these polynomials can then be used to evaluate the needed inner products efficiently. The 
term A> is evaluated as follows, where for simplicity we assume that the terms are 

functions of just a single variable, r. 

(1) Evaluate ^ at the collocation points, r,-, t = 1, . . . ,iV. 

(2) Evaluate ^ there also. 

(The derivatives of are evaluated spectrally.) 

(3) Evaluate ^ at the collocation points. 

(4) Express each component of ^ as a series of the appropriate polynomials. 

(5) Evaluate the integrals <fti, A> utilizing the orthogonality properties of these polynomials. 

When using Tchebychev polynomials to form the expansion and projection vectors, the 
transformations needed in step (4) can be very efficiently carried out (in 0{NlogN) operations) 

by means of Fast Fourier Transforms (see, for example, (Canute et al.)); however, it is then 

1 

important to include the weight function, (1 - r^) ^ (assuming r e [-1,1]), in the definition of 
the projection vectors so that we can carry out step (5) efficiently (Moser &; Moin [1984]). We 
also mention in passing that to avoid aliasing errors it is necessary to use more collocation points 
than the number of modes, N, in the expansion. 

For the time discretization scheme, the viscous terms are generally advanced by means of 
the trapezoidal rule (Crank-Nicholson) with an Adams-Bashforth scheme being used for the non- 
linear term, /. Such a scheme is indicated in the following equation where the superscript denotes 
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the time level and St is the step size being used: 

= . 1 ( 3 ^ - . (4.1) 

The velocity is represented by a truncated expansion, with the coefficents now depending on 

the time, t. Projections are taken with vectors to get the following set of linear equations to 
solve for the coefficents, a at the new time level, j+1: 

I “ I f ( 

where C and A are essentially the same matrices we had before arising respectively from the 
viscous terms and the non-orthogonality of the expansion and projection vectors terms while F 
represents the projection of /, 

F = </, A>. (4-3) 

Moser & Moin [1984] used a similar scheme to compute several states of the infinite Taylor- 
Couette flow successfully and went on to simulate turbulent flow in a mildly curved channel. For 
example, using just nine Fourier modes in the z direction and eleven Tchebychev polynomials in 
the r direction they found that the critical Reynolds number for transition to Taylor vortices was 
at Re = 185 for cylinders with an inner to outer ratio, r} = 0.95 and at Re = 68.2 when rj = 0.5. 
These figures are in excellent agreement with the analysis of DiPrima & Eagles [1977] for both 
the narrow and wide gap problems. They also carried out a series of three dimensional Taylor- 
Couette flow calculations and had excellent agreement with experimental results for Reynolds 
numbers up to 1300 with rj = 0.877 (for modulated wavy vortex flow) using thirty one Fourier 
modes in the z direction, fifteen in the 6 direction, and thirty three Tchebychev polynomials in 


the radial directon. 
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5. Conclusions and Outlook 


Divergence free expansions eliminate many of the problems that a numericist must face 
when discretizing the equations governing the flow of incompressible fluids. Use of such an 
expansion ensures that the computed velocity will be exactly solenoidal and also enforces all the 
necessary boundary conditions on that vector, including the no-slip condition on the solid boun- 
daries. Projection with inviscid vectors eliminates the yp term from the governing equation and 
we can recover the pressure, if it is needed, once an adequately resolved velocity field has been 
obtained. Finally , as each component of the velocity is expressed as a sum of polynomials which 

are solutions to a singular Sturm Liouville problem, we can expect exponential convergence for a 
smooth u. 

There are some problems with this method also. In common with all global expansion 
methods It IS not possible to consider complex physical domains directly - these must be mapped 
onto some simpler configuration, which may not always be possible, or else split into simple sub- 
domains, leaving us with the problem of connecting the solutions across the interfaces. This 
problem is by no means specific to the divergence free expansion technique but rather is inherent 
in all spectral methods. It must also be admitted that implementation of the vector expansions 
in a computer program is far from easy, requiring as it does, a significant amount of algebraic 
manipulation to subdivide the matrices that are used into submatrices that are independent of 
the parameters occurring in the problem. A symbolic language such as MACSYMA can be used 
to somewhat alleviate these latter difficulties. 

Current research that utilizes divergence free expansions includes a simulation of the finite 
Taylor-Couette system with the aid of vectors that are aperiodic in both the radial and axial 
directions. We see this method as an excellent tool for the direct simulation of incompressible 
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flows, and it should have many applications in investigations of the fundamental problems of 
incompressible fluid mechanics. 
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